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ABSTRACT 

We study the formation of globular clusters in a Milky Way-size galaxy using a high-resolution cosmo- 
logical simulation. The clusters in our model form in the strongly baryon-dominated cores of supergiant 
molecular clouds in the gaseous disks of high-rcdshift galaxies. The properties of clusters are estimated 
using a physically-motivated subgrid model of the isothermal cloud collapse. The first clusters in the 
simulation form at z « 12, while we conjecture that the best conditions for globular cluster formation 
appear to be at z ~ 3 — 5. Most clusters form in the progenitor galaxies of the virial mass Mh > 10 9 M 
and the total mass of the cluster population is strongly correlated with the mass of its host galaxy: 
Mqc = 3 x 10 6 M©(M h /10 n Mo) 11 . This corresponds to a fraction - 2 x KT 4 of the galactic baryons 
being in the form of globular clusters. In addition, the mass of the globular cluster population and the 
maximum cluster mass in a given region strongly correlate with the local average star formation rate. 
We find that the mass, size, and metallicity distributions of the globular cluster population identified 
in the simulation are remarkably similar to the corresponding distributions of the Milky Way globulars. 
We find no clear mass-metallicity or age-metallicity correlations for the old clusters. The zero-age mass 
function of globular clusters can be approximated by a power-law dN/dM oc M~ Q with a 2, in agree- 
ment with the mass function of young stellar clusters in starbursting galaxies. We discuss in detail the 
origin and universality of the globular cluster mass function. Our results indicate that globular clusters 
with properties similar to those of observed clusters can form naturally within dense gaseous disks at 
z > 3 in the concordance ACDM cosmology. 

Subject headings: cosmology: theory-galaxies: formation-globular clusters: formation - methods: 
numerical 



1. INTRODUCTION 

More than seventy years ago in a monograph entitled 
Star Clusters Harlow Shapley (1930) wrote: "It is encour- 
aging to see how fragile and futile are the majority of as- 
tronomical theories and speculations... for the futility of 
speculations emphasizes the importance and durability of 
observations and indicates the steady progress of the sci- 
ence." These words are particularly true for the models of 
globular cluster (GC) formation. Extensive observational 
surveys of globular cluster systems in the Milky Way and 
other galaxies have been compiled during the past two 
decades (e.g., Harris 2001). At the same time, despite a 
wide variety of proposed models, a self-consistent scenario 
of globular cluster formation is yet to be constructed. 

The existing models can be classified into four broad 
categories. In the primary models globular clusters are 
envisioned to have formed soon after recombination, with 
masses determined by the cosmological Jeans mass (Pee- 
bles & Dicke 1968; Peebles 1984). In the secondary mod- 
els, globular clusters are assumed to appear during the 
early stages of galaxy formation. For instance, Fall & Rees 
(1985) pointed out that thermal instability in hot gaseous 
halos of young galaxies can naturally lead to the conden- 
sation of globular cluster-size clouds. Several other trigger 
mechanisms operating during galaxy formation, such as 
the shock compression and collisions of primordial molec- 



ular clouds, were also explored (Gunn 1980; Burkert et al. 
1992; Murray & Lin 1992; Larson 1996; Harris & Pudritz 
1994; Cen 2001). 

Models of the third class correspond to the relatively 
recent stages of galaxy formation. Schweizer (1987) and 
Ashman & Zepf (1992), for example, proposed a model 
of GC formation in the gas-rich mergers of disk galax- 
ies. There mergers perturb, compress, and shock the in- 
terstellar medium which creates the very high pressure re- 
gions conducive to GC formation. Accordingly, this model 
predicted that young clusters should be present in merg- 
ing and interacting galaxies. These predictions were suc- 
cessfully confirmed by HST observations (Whitmore & 
Schweizer 1995; Holtzman et al. 1996; Whitmore et al. 
1999; Zepf et al. 1999). 

The division between the second and third classes of 
models is somewhat blurred. In the current hierarchi- 
cal structure formation paradigm, galaxy formation is a 
continuous process of merging and accretion. The fourth, 
most recent class of models thus incorporates the elements 
of all previous classes within the hierarchical framework 
(Cote et al. 2000, 2002; Beasley et al. 2002; Gnedin 2003). 
Globular clusters in these models are assumed to form 
in young disks which undergo frequent mergers. At this 
point, the models are largely phcnomcnological and char- 
acterize the cluster formation using somewhat ad hoc recipes, 
limiting their predictive power. Nevertheless, they have 
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been fairly successful in reproducing many properties of 
the observed GC populations. The comparisons of various 
models to observations can be found in Harris (2001) and 
Gnedin et al. (2001). 

The main obstacle to building a realistic and self-consistent 
physical model of globular cluster formation has always 
been the uncertainty in the initial conditions. In fact, 
all of the models mentioned above would produce star 
clusters in environments where the conditions agree with 
the model assumptions. Another unknown is the connec- 
tion between globular clusters and galaxies. On the one 
hand, the largest globular clusters have masses compara- 
ble to those of dwarf spheroidal galaxies (~ 10 7 M Q ). On 
the other, globular clusters do not seem to have extended 
dark matter halos (e.g., Moore 1996) and in that respect 
differ fundamentally from galaxies. There is also signifi- 
cant disparity between the densities, velocity dispersions, 
and structural parameters of dwarf galaxies and globular 
clusters (Kormendy 1985). In order to understand these 
differences we need a self-consistent model which ties the 
formation of globular clusters to the realistic formation 
and evolution of their parent galaxies. 

The theory of hierarchical galaxy formation has matured 
in the last decade, motivated by the theory of inflation, 
guided by observations, and aided by elaborate numer- 
ical simulations. In hierarchical models, galaxies form 
via gravitational instability from small-amplitude initial 
Gaussian fluctuations with well-defined statistical proper- 
ties. Recently, this scenario has been spectacularly con- 
firmed by the CMB anisotropy measurements and other 
cosmological probes (e.g., Spergel et al. 2003). The spa- 
tially flat cosmological model dominated by the dark en- 
ergy and cold dark matter (ACDM) favored by observa- 
tions provides a solid framework for the theory of globular 
cluster formation. 

Cosmological simulations follow the hierarchical build- 
up of galaxies self-consistently starting from the well-defined 
initial conditions. These simulations are now reaching 
the level of sophistication and dynamic range sufficient 
to study the formation and dynamics of giant molecular 
clouds in galactic disks. Therefore, we can address the 
formation of the proto-cluster clouds without resorting to 
phenomenological parameterization and directly study the 
details of when, where, and how globular clusters formed. 

The main goal of this work is to study the formation of 
globular cluster populations in the hierarchical scenario us- 
ing a very high-resolution cosmological simulation. Based 
on observational evidence, we assume that clusters form in 
dense isothermal cores of the super giant molecular clouds 
ubiquitous in high-rcdshift galactic disks. In addition, we 
use a simple model of isothermal collapse to derive the 
properties of stellar clusters that would form in such cores. 
We then compare the derived properties of model globular 
clusters with those of the GC population in the Galaxy 
as well as with the populations of young GCs in external 
galaxies. 

Many decades ago, Harlow Shapley used the distribution 
of globular clusters to greatly expand and re-define the 
structure of our Galaxy. It is only fitting that we now 
apply our understanding of galaxy formation to predict 
and explain the properties of globular clusters. 



2. NUMERICAL SIMULATIONS 

2.1. Numerical techniques and Physical Processes 

The simulation presented in this paper was performed 
using the Eulerian gasdynamics+A-body Adaptive Re- 
finement Tree (ART) code. This code is based on the cell- 
based approach to adaptive mesh refinement developed by 
Khokhlov (1998). The algorithm uses a combination of 
multi-level particle-mesh (Kravtsov et al. 1997; Kravtsov 
1999) and shock-capturing Eulerian methods (van Leer 
1979; Colella & Glaz 1985) to follow the evolution of dark 
matter (DM) and gas, respectively. High dynamic range 
is achieved by applying adaptive mesh refinement both to 
the gasdynamics and gravity calculations. 

Several physical processes critical to various aspects of 
galaxy formation are implemented in this code: star for- 
mation; metal enrichment and thermal feedback due to the 
supernovae type II and type la (SNII/Ia); self-consistent 
advection of metals; metallicity- and density-dependent 
cooling and UV heating due to the cosmological ionizing 
background, using the cooling and heating rates tabulated 
in the temperature range 10 2 < T < 10 9 K for a grid of 
densities, metallicities, and UV intensities using Cloudy 
(ver. 96b4, Ferland et al. 1998). In the present simula- 
tions we set a minimum temperature of T m ; n = 300 K. 
The cooling and heating rates take into account Comp- 
ton heating/cooling of plasma, UV heating, atomic and 
molecular cooling. While the detailed implementation of 
these processes is described elsewhere (Kravtsov 2003a, b), 
below we summarize the details crucial to this study. 

We use a "constant efficiency" star formation prescrip- 
tion. Namely, the stars are formed with a constant timescale 
t* so that the star formation rate is proportional to the lo- 
cal gas density, oc p g . This prescription is motivated by 
observations of the star forming regions (e.g., Young et al. 
1996; Wong & Blitz 2002) and appears to reproduce the 
Schmidt-like law of star formation on kpc scales (Kravtsov 
2003b). Star formation is allowed to take place only in the 
coldest and densest regions, T < T$f and p g > psf, but 
no other criteria (like the collapse condition V • v < 0) 
are imposed. We use r* = 4 Gyr, Tsf = 9000 K, and 
Psf = 1-64 M Q pc~ 3 or the atomic hydrogen number den- 
sity of nn = 50 cm -3 . The adopted values of Tsf and psf 
are quite different from the typical temperatures and den- 
sities of star forming molecular cores: T < 30 — 50 K and 
7iH > 10 4 cm~ 3 . They are, however, more appropriate for 
the identification of star forming regions on ~ 100 pc scales 
which are resolved in the present simulation. In practice, 
Tsf is not relevant because most of the gas with p > psf 
is at temperatures of just a few hundred degrees Kelvin. 

Each newly formed stellar particle is treated as a single- 
age stellar population and its feedback on the surround- 
ing gas is implemented accordingly. The word feedback 
is used here in a broad sense to include the injection of 
energy and heavy elements (metals) via stellar winds and 
supernovae, and the secular stellar mass loss. Specifically, 
we assume that the stellar initial mass function (IMF) is 
described by the Miller & Scalo (1979) functional form 
with stellar masses in the range 0.1 — 100 M Q . All stars 
with m* > 8 Mq deposit 2 x 10 51 ergs of thermal en- 
ergy and a mass /zrn* of heavy elements in their par- 
ent cell (no delay of cooling is introduced in these cells). 
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The metal fraction is fz = min(0.2, 0.01m* — 0.06), which 
crudely approximates the results of Woosley & Weaver 
(1995). In addition, the stellar particles return a fraction 
of their mass and metals to the surrounding gas at a sec- 
ular rate mi oss = rn„ Co(t — tbirth + To) with Co = 0.05 
and T = 5 Myr (Jungwiert et al. 2001). The released 
metals are advected along with the gas. The code also 
accounts for SNIa feedback assuming a rate that slowly 
increases with time and broadly peaks at the population 
age of 1 Gyr. We assume that the fraction 5 x 10~ 3 of mass 
in stars between 3 and 8 M Q explodes as SNIa over the 
entire population history and each SNIa dumps 2 x 10 51 
ergs of thermal energy and ejects 1.3 M Q of metals into 
parent cell. For the assumed IMF, 75 SNII (instantly) 
and 11 SNIa (over several billion years) are produced by a 
10 4 M© stellar population. 

2.2. Simulation Parameters 

The simulation we use in our analysis follows the early 
( z ^ 3) stages of evolution for a galaxy of typical mass: 
« 10 12 ft, _1 M Q at z = 0. At the analyzed epochs, the 
galaxy has already built up a significant portion of its final 
mass: 1.3 x lO 10 ^ 1 M Q at z = 9 and 2 x lO 11 /!" 1 M Q at 
z = 4. The total galaxy mass, Mh, is defined as the mass 
enclosed within the radius of the average density equal 
340 times the mean matter density. The simulation starts 
from a random realization of the Gaussian density field 
at z = 50 in a periodic box of 6/i _1 comoving Mpc with 
the power spectrum (Hu & Sugiyama 1996) appropriate 
to the flat ACDM model: fl = 1- Ct A = 0.3, Q, b = 0.043, 
h = _ff /100 = 0.7, n s = 1, and <7g = 0.9. The parameters 
have their usual meaning and are consistent with recent 
cosmological constraints (e.g., Spergel et al. 2003). 

To increase mass resolution, a low resolution simulation 
was run first and a galactic-mass halo was selected. A 
lagrangian region corresponding to five virial radii of the 
object at z = was then identified at z — 50 and re- 
sampled with additional small-scale waves (Klypin et al. 
2001). The total number of DM particles in the high- 
resolution lagrangian region is 2.64 x 10 6 and each particle 
mass is muM = 9.18 x 10 5 ft- _1 M Q . Outside the high- 
resolution region the matter distribution was sampled with 
« 3 x 10 5 higher mass particles. 

As the matter distribution evolves, the code adaptively 
and recursively refines the mesh in the high density re- 
gions. Initially, a uniform 64 3 grid covered the entire com- 
putational box. The lagrangian region, however, was al- 
ways unconditionally refined to the third refinement level, 
corresponding to an effective grid size of 512 3 . Beyond the 
third level, a mesh cell was tagged for refinement if its gas 
or DM mass exceeded 0.125 and 0.0625 times the mean 
mass expected for the average density in each component 
in the zeroth level (i.e., uniform grid) cell, respectively. 
The refinement follows the collapse of 1.2 x 10 6 /i _1 M 
(gas) and 3.7 x 10 6 /i _1 M (DM) mass elements in a quasi- 
lagrangian fashion. These masses can be loosely consid- 
ered as gas mass resolution until the maximum level is 
reached beyond which refinements arc not done. In the 
run we use this level is set to Z max = 9 and is reached by 
z w 10. Beyond this, the notion of mass resolution for gas 
is not well defined because gas is represented as a contin- 
uous medium on an Eulcrian mesh. Once the maximum 



refinement level is reached, the mass per cell then is no 
longer constant but reflects the local gas density. For ex- 
ample, cells of the 9th level have gas densities spanning 
the range of more than six orders of magnitudeas the in- 
terstellar medium is multi-phase with tenuous hot gas and 
very dense cold gas occupying different regions (see, e.g., 
Fig. 9 below). 

The spatial resolution of the simulation is thus time- 
dependent. As the density increases, additional refiniment 
levels are added to keep the mass per cell approximately 
constant. The maximum allowed refinement level Z max was 
set to nine and this level was reached at z « 10. In 
the simulation we present, the physical size of the max- 
imum refinement cell is ~ 28, 20, 26, 37, and 45/i _1 pc 
at z — 12, 9, 6, 4, and 3, respectively. Thus, the change 
over the analyzed range of epochs is not very large. A 
total of s=a 1.1 x 10 7 mesh cells was used at z — 4 with 
« 2.5 x 10 5 of them at refinement levels of 8 and 9. The 
high-density cold star forming disks within DM halos were 
refined to l max — 9. The physical size of mesh cells was 
Ax t = 26.16 [10/(1 + z)}2 9 - 1 pc, where I is the cell's level 
of refinement. Each refinement level was integrated with 
its own time step Ati = At 2~ l w 2 9 ~ z x 2 x 10 4 yr, where 
Ato < 10 7 yr is the global time step on the zeroth level 
set using the Courant-Friedrichs-Levy condition. 

2.3. Identification of the Globular Cluster Formation 

Sites 

Although the resolution achieved in the simulation is 
very high in the disk region, it is still insufficient to resolve 
the formation of stellar clusters. The resolution, however, 
is sufficient to identify the potential sites for GC forma- 
tion. The cores of giant molecular clouds in high-redshift 
galaxies are the natural candidates (Harris & Pudritz 1994; 
McLaughlin & Pudritz 1996) for such sites. Numerical 
simulations of Nakasato et al. (2000) show that globular 
clusters with realistic masses and sizes can indeed form in 
such cores. We therefore adopt this picture and identify 
the cores of dense gaseous clouds in simulated disks as the 
sites of globular cluster formation. 

We identify giant molecular clouds using the following 
algorithm. All mesh cells with gas densities greater than 
a certain threshold density, p mc , are selected and sorted 
into a list of increasing density. The highest density cell 
and all of its immediate neighbors are then included in 
the first cloud. This cell is labeled as the core cell of the 
cloud. The next highest density cell in the list is then 
considered. If it happens to be already included in the 
first cloud, all of its immediate neighbor cells are then also 
included in the first cloud. If, on the other hand, the cell 
is not part of a cloud yet, it is assigned to a new cloud 
and is labeled as its core. The procedure repeats until the 
list of cells is exhausted. The algorithm is thus somewhat 
similar, but not equivalent, to the well-known friends-of- 
friends clustering algorithm. The current implementation 
will break a region into separate clouds for each density 
peak rather than combining several peaks into the same 
cloud. 

We explored several values of the threshold density p mc = 
1 — 50 M Q pc~ 3 . Although the cloud masses grow and ever 
smaller clouds are included in the catalog as the threshold 
is decreased, in the mass range relevant for the GC iden- 
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Fig. 1. — Gas density (left panel) and temperature (right panel) in the most massive disk at z = 4. The projected density is in cm -2 , while 
the density-weighted average temperature is in degrees Kelvin. The vectors in the right panel show gas velocities; the thick vertical vector in 
the lower left corner of the panel corresponds to 200 kms" 1 . The density and temperature are projected over a 3.5 kpc slice centered on the 
cell of the maximum gas density (the center of the plot). The figure shows a nearly face-on disk with prominent spiral arms in the process 
of very active accretion and merging. In our model the globular clusters form in the densest regions of the disk corresponding to the darkest 
knots in the temperature map. The globular clusters identified in the disk at this epoch are shown by circles in the left panel. The radius of 
the circles corresponds to the mass of each cluster. 



tification the same cores are recovered for all p mc - For our 
analysis below we choose the cloud catalog with the fidu- 
cial value p mc — 1M Q pc~ 3 . This corresponds to the gas 
number density « 40 cm -3 and pressure > 10 4 K cm~ 3 . 
Note that at these densities the gas temperature is at the 
lowest value allowed in the simulation: T m ; n = 300 K. 

2.4. The Subgrid Model 

In order to derive the properties of star clusters forming 
in molecular cloud cores, we complement the simulations 
with a physical description of the gas distribution on a 
subgrid level. Observations of star forming regions in the 
solar neighborhood (e.g., Elmegreen 2002) show that clus- 
tered star formation proceeds inevitably when the cores 
of molecular clouds reach a critical density > 10 5 cm~ 3 . 
Formation of young massive clusters is also accompanied 
by high external pressure, P > 10 7 K cm -3 . Within 
such cores, star clusters form with a locally-high efficiency 
(Geyer & Burkert 2001; Kroupa et al. 2001; Kroupa & 
Boily 2002): 



which is required in order to produce gravitationally-bound 
clusters. On the theoretical side, analytical models and 
numerical simulations also indicate that dense bound clus- 
ters can form quickly and efficiently in the cores of giant 
molecular clouds of moderate metallicity (e.g., Harris & 
Pudritz 1994; McLaughlin & Pudritz 1996; Nakasato et al. 
2000). In the models of molecular clouds including thermal 
support, turbulence, and magnetic fields star formation 
proceeds rapidly, on one or two dynamical times (Pringle 



1989; Ostriker et al. 2001; Bate et al. 2003). Based on 
these observational and theoretical results, we implement 
the following phenomcnological subgrid model. 

At high densities in the central cells of the identified gas 
clouds the cooling time is much shorter than the dynam- 
ical time, and therefore the cells must be isothermal. In- 
deed, the resolved structure of the molecular clouds in the 
simulation has an isothermal profile, p g oc r~ 2 . We thus 
extrapolate this profile inside the central unresolved cell. 
(Such isothermal structure is also predicted by the simula- 
tions of the collapse of cloud cores by Nakasato et al. 2000). 
What we measure in the simulation is the cell-averaged gas 
density, p cc \\ = p av (< R cc n), where i? C eii is the cell radius 
(a half of the cell dimension) . For an isothermal profile the 
average density within i? ccl i is p av (< R cc ii) = 3p g (R cc n). 
We can thus derive the inner density profile as 

p g {r < R cell ) = i p coU {^-^j ■ ( 2 ) 

We assume that a single cluster forms within the cloud 
core on a dynamical (free-fall) time at the densities higher 
than the critical, p cs f. To distinguish from the field star 
formation denoted by the subscript 'SF' in §2.1, here the 
subscript 'csf stands for 'clustered star formation'. By 
choosing a density threshold we postulate that only a high- 
density tail of the gas distribution produces compact mas- 
sive clusters, whereas the rest (and most) of the gas par- 
ticipates in the formation of field stars or open clusters. 
This scenario is quite natural. High gas densities are re- 
quired to match the observed stellar densities, which are 
highest in globular clusters. A single value of the density 
threshold, of course, is not sufficient to describe complex 



physics of molecular clouds but it can be very useful in 
denning the properties of the clusters, as we demonstrate 
below. 

The radius of the molecular core going into clustered 
star formation i? cs f is determined by the condition p g (R cs f) - 
p cs f. All gas within i? cs f is converted into stars with the 
efficiency e: 

M = e M g {R csl ) = e 4irp csf i? 3 sf . (3) 

The fraction 1 — e of the core mass remaining in the gas 
phase will be expelled from that region, following the for- 
mation of UV-bright O and B stars. As a result of the 
gradual loss of the remaining gas, the star cluster expands 
almost adiabatically (e.g., Geyer & Burkert 2001; Boily & 
Kroupa 2003) such that its final size is 
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The resulting average density of the clusters is 
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We use the fiducial values e = 0.6 and p cs f = 10 M Q pc~ 3 , 
which gives p « 4 x 10 3 M© pc~ 3 . It is close to the median 
density at the half-mass radius for the Milky Way globular 
clusters, which is 3 x 10 3 M pc -3 . 

Although the expressions for the mass and size of star 
clusters are linked to the resolution-dependent cell proper- 
ties (pceii, -Rceii), the parameters of individual clusters are 
almost insensitive to changes in the resolution. This is be- 
cause the cores of dense molecular clouds in which globular 
clusters form are isothermal. Thus, when the cell size i? C oii 
changes, the cell density adjusts as p cc \\ oc R~^, leaving 
the radius R cs { and the cluster radius R and mass M un- 
changed (see eqs. [3] and [4]). This is indeed what we find 
in the test presented below in § 5. When we repeat the 
analysis decreasing the level of refinement (Z max — 1), the 
masses of individual clusters do not change significantly. 

In the following, we present the properties of the model 
clusters based on the density criterion. For completeness, 
in §6.3 we discuss and evaluate alternative subgrid models 
and show that our model works best in reproducing the 
observed properties of clusters. 

Note that the mass function of globular clusters at birth 
will be significantly modified by the effects of dynamical 
evolution. As we discuss in §6.1, low-mass and low-density 
clusters are preferentially dissolved by the combined effects 
of two-body relaxation, tidal shocking, dynamical friction, 
and stellar evolution. High-mass clusters (M > 10 5 M©), 
on the other hand, preserve their mass function and trace 
the initial distribution. 

3. RESULTS 

We analyze the simulation outputs at twelve epochs be- 
tween redshifts z — 11.8 and 3.35, identifying the cores 
of the giant molecular clouds and computing properties of 
the model globular clusters, as described in § 2.3 and 2.4. 
The time intervals between the outputs are in the range 
~ 100 — 300 Myr. Due to limited computational resources, 
the simulation was run only until z — 3.3. 

3.1. Spatial distribution of globular clusters at high 
redshifts 






Fig. 2. — The identified globular clusters within the global dis- 
tribution of dark matter at z = 7 and 2 = 4. The view is centered 
on the largest galaxy in the simulation and shows lh~ 1 Mpc region 
(comoving). The gray-scale colored particles represent the dark mat- 
ter, while white circles in the centers of some halos show locations of 
the globular clusters identified in the simulation. Note that massive 
halos contain multiple clusters in their centers. The dark matter 
particles are colored according to the local density on a logarithmic 
stretch. 



Before discussing the detailed properties of model clus- 
ters, we first consider their spatial distribution. Figure 1 
shows the density and temperature maps, as well as the 
velocity field of gas, in the region of the most massive 
disk at z = 4. The cooling of the gas in these regions is 
very fast and the cold gas always settles in a thin disk. 
Frequent interactions drive strong spiral density waves in 
the gaseous disks, which fragment into separate molecular 
clouds. Globular clusters identified within those clouds 
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are represented by circles in the density map. The figure 
shows that clusters in our model form in the dense cold re- 
gions, generally tracing the spiral arms of the galaxy. The 
morphology of the distribution is very similar to that of 
young stellar clusters observed in starburst galaxies (e.g., 
Whitmore & Schweizer 1995; Zhang et al. 2001). 

At this epoch the parent halo of the disk experiences 
frequent mergers and vigorous accretion of fresh gas. The 
two cold knots outside the spiral arms, for instance, are 
the small-mass satellite galaxies in the process of merging 
with the central halo. The dense gas in these satellites 
could have been compressed by the external pressure and 
shocks accompanying their collision with the disk. The 
galaxy as a whole exhibits frequent bursts of star formation 
associated with minor and major mergers. These mergers 
destroy the short-lived disks and scatter away young stars 
and star clusters in a spheroidal halo. The cold gas, on 
the other hand, always falls back to form a new thin disk. 

Figure 2 shows the spatial distribution of model globular 
clusters within the large-scale structure formed in our sim- 
ulation at z — 7 and z = 4. The distribution of dark mat- 
ter is typical of hierarchical models. Visually, it is domi- 
nated by prominent filaments on large scales and hundreds 
of dense dark matter halos tracing these filaments on small 
scales. The figure shows that parent halos of globular clus- 
ters are tracing the skeleton of the large-scale structure. 
They concentrate in the densest regions of the filaments 
close to the central massive object. In other words, the dis- 
tribution of halos containing clusters is highly biased with 
respect to the overall distribution of matter. This bias 
is especially pronounced at z = 7. The highly clustered 
distribution of halos at the early epochs is a generic fea- 
ture of the hierarchical models, in which objects of galactic 
mass correspond to the relatively high peaks in the initial 
Gaussian density field. This property can be extremely 
important for explaining the present distribution of glob- 
ular clusters in the halo of the Milky Way. The high spatial 
bias of globular clusters at early epochs would result in the 
more concentrated radial distribution of globular clusters 
compared to the dark matter today (West 1993) and in the 
preferential location of higher-metallicity clusters towards 
to the center, in agreement with observations (Djorgovski 
& Meylan 1994; van den Bergh 2003). 

Note that although the high-rcdshift globular clusters 
form in dense gaseous disks, the subsequent accretion of 
their parent galaxies along filaments will lead to tidal strip- 
ping and disruption. For example, analysis of the evolu- 
tion of Milky Way size progenitors (Kravtsov, Gnedin, & 
Klypin 2004) shows that most of the dwarf-size systems 
that are located within < 3 virial radii from the progenitor 
at z > 4 accrete early and are disrupted before present day 
epoch. This includes most of the objects hosting globular 
clusters in Figure 2. The disrupted systems form diffuse 
dark matter halo and contribute to the stellar halo of the 
host (Bullock, Kravtsov, & Weinberg 2001). Their clus- 
ters would share the fate of the stripped stars and should 
therefore have spatial distribution at z — similar to the 
stellar halo stars. Direct observational evidence of disrup- 
tion is provided by the extended tidal tails around glob- 
ular clusters (e.g., Odenkirchen ct al. 2003; Sohn et al. 
2003), around dwarf satellite galaxies (Freeman & Bland- 
Hawthorn 2002, and references therein), and the possible 
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Fig. 3. — Accretion history in the central cell of the main progen- 
itor halo. Upper panel: Gas number density. The densities of dark 
matter and stars in this region are an order of magnitude smaller. 
Lower panel: Iron abundance of the gas with respect to the solar 
value. The metallicity is due to SNII enrichment, as the contribution 
of SNIa is negligible at these epochs. 

association between the two (Lynden-Bell & Lynden-Bell 
1995). We will investigate the dynamical evolution and 
the present-day spatial distribution of the GC population 
in our model in a future study. 

3.2. Molecular clouds in the dense high-redshift disks 

The globular clusters in our model form in the high- 
density cores of giant molecular clouds of the high-rcdshift 
galaxies (Fig. 1). It is therefore important to consider the 
properties of the molecular clouds in connection to the 
properties of globular clusters and host galaxies. Figure 
3 shows the evolution of density and metallicity in the 
central cell of the most massive disk shown in Figure 1. 
Although this cell has the highest gas density and is lo- 
cated at the bottom of potential well, the overall evolution 
is common for all cells. The gas density exhibits several 
prominent peaks associated with the fast episodes of ac- 
cretion. Cold metal-poor gas is delivered to the center of 
the disk both by merging of smaller galaxies and by direct 
accretion of gas along a filament that reaches inside the 
disk region. The rapid increase of the density and pres- 
sure in the cell during the accretion events can trigger the 
collapse of the molecular cloud. The gas density saturates 
at lower redshifts (z < 5) as the accretion on the center of 
the disk slows down. 

The bottom panel of Figure 3 shows the evolution of 
metallicity due to the SN ejecta (the contribution of SNIa 
to the metallicity is negligible at these epochs) . The metal- 
licity quickly increases to about 10% of solar and then 
evolves slowly. Note that the events of accretion of the 
fresh low-metallicity gas may lower the mean metallicity, 
even in the very central region. If a series of globular clus- 
ters forms between z — 8 and 4 in this region, the younger 
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ones are not necessarily more metal-rich than the older 
ones. The age difference of these clusters would however 
be less than 2 Gyr. 

Overall, the galaxies in the simulation exhibit a well- 
defined correlation between the stellar mass and the av- 
erage metallicity of stars, Z oc M°- 5 , similar to the cor- 
relation observed in nearby dwarf galaxies (Dekel & Woo 
2003). There is also a significant spread in gas metal- 
licity even within a single object, which indicates that 
mixing of metals is rather inefficient. The wide range of 
gas metallicity in star forming regions eliminates any clear 
age-metallicity correlation for high-redshift clusters. For 
instance, stars formed at the same epoch can have metal- 
licities different by up to two orders of magnitude. This 
may also at least partially explain the well-known "second 
parameter problem" (e.g., Carney 2001). 

The efficiency of GC formation, i.e. the ratio of the glob- 
ular cluster mass to the mass of the molecular, baryonic, 
or dark matter system containing it, depends on the aver- 
aging scale. Within the cores of giant molecular clouds the 
local efficiency is of the order unity (§ 2.4). Averaged over 
the whole molecular cloud though, the efficiency is much 
lower because most of the molecular gas is not participat- 
ing in star formation at any given time. When compared 
with the total gas and/or dark matter mass in the galaxy, 
the efficiency decreases by another order of magnitude. 
There are thus various types of globular cluster formation 
efficiencies, which we consider in turn. 

The detailed properties of the simulated molecular clouds 
depend on the threshold density, p mc , used to define the 
cloud boundary (see § 2.3). This boundary can be thought 
of as an external tidal limitation. The mass and size of the 
cloud increase with the decreasing threshold density. The 
cloud-scale efficiency of globular cluster formation, which 
we define as ecc = M/M mc , varies accordingly. We find 
that the average efficiency is about 10~ 2 for p mc = 50 M 
pc~ 3 , is in the range 10~ 3 - 10~ 2 for p mc = 10 M pc~ 3 , 
and 10~ 4 - 10~ 3 for p mc = 1 M pc~ 3 . The estimated 
masses of globular clusters, on the other hand, depend 
only the properties of the cloud cores and are insensitive 
to the changes in the external boundary condition. 

For the massive globular clusters with M > 3 x 10 5 M 
and the associated massive molecular clouds in our fiducial 
model with p mc = 1 M pc~ 3 , the formation efficiency is 
roughly constant: ecc ~ 10~ 3 . However, if we include 
lower mass clusters we find an anti-correlation with the 
cloud mass (Spearman correlation coefficient r s = —0.35). 
In the range 10 5 M < M mc < 1O 8 M , the relation is 
lge GC = -1.6 - (0.22 ± 0.02) lg(M mc /M ). Overall, the 
numerical values of the efficiencies we obtain are in good 
agreement with observations (Harris & Pudritz 1994). 

3.3. The global efficiency of globular cluster formation 

An important measure of the efficiency of globular clus- 
ter formation is the total mass of clusters, Mqc, within 
a parent galactic halo. For example, in giant elliptical 
galaxies the ratio of the total cluster mass to the mass 
of stars plus the hot X-ray emitting gas is roughly con- 
stant, £^ c = Mqc /Mb « 0.0026 ± 0.0005 (McLaugh- 
lin 1999). This parameter can be thought of as the ef- 
ficiency of the conversion of baryons into globular clus- 
ters. In massive objects the baryon mass Mb relates to 
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Fig. 4. — The mass of the globular cluster system within a given 
halo vs. the total mass of its parent halo, combined over all analyzed 
epochs. Dots show the average in bins of width A lg = 0.2, while 
the error bars show the ltr deviations within the bin (not the error 
of the mean). The solid line is the least-squares fit with the slope 
d\gM GC /dlgM h = 1.13 ± 0.08. 



the total galaxy mass Mh via the universal baryon frac- 
tion, Mb/Mh ~ fb ~ 0.14. Thus perhaps even more fun- 
damental is the ratio of the globular cluster mass to the 
total galaxy mass: £q C = Mqq/M^. 

Figure 4 shows the sum of the globular cluster masses 
in each halo versus the progenitor galaxy mass at the time 
of GC formation. There is a well-defined correlation of the 
form 

/ M \ 1.13±0.08 

MGC = 3 - 2xl ° 6M Hl0n^J ® 

albeit with scatter. The global efficiency £q C is therefore 
only weakly dependent on the galaxy mass. For most halos 
harboring massive clusters we find £q C — (2 — 5) x 10~ 5 . 
The global baryon efficiency is in the range £q C = (2 — 
3) x 10~ 4 , and it scales with the galaxy mass as 

0.25±0.12 



M GC 



2.5 x 10" 



M h 



(7) 



M b V 1 O 11 M 0/ 

While these values are lower than those found by McLaugh- 
lin (1999), they are in fact appropriate for a spiral galaxy 
like our own. McLaughlin (1999) derived £q C = 0.0027 
for the Galaxy taking into account only the mass of the 
stellar spheroid. We are interested here in scaling with the 
total mass and therefore take the most recent estimate of 
the virial mass of the Milky Way M h « 10 12 M (Klypin 
et al. 2002). The mass of the observed globular clusters 
is Mqc = 5.2 x 1O 7 M (Harris 1996), and as McLaugh- 
lin (1999) argues, this mass cannot differ from the initial 
mass by more than 25%. Thus the global efficiency for 
the Galaxy is £q C > 5 x 10~ 5 , and the baryon efficiency 
£q C ~ £Q C //b ~ 4 x 10~ 4 . Both of these estimates agree 
with our derived correlations, eqs. (6) and (7). 
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Fig. 5. — The build-up of the initial mass function of globular 
clusters. Dotted, dashed, long-dashed, and solid histograms show 
cumulative distributions at z = 10, 7, 4, and 3.3, respectively. The 
straight dashed line shows a power-law, N oc M~ a , with the slope 
a = 2. Note that this is the mass function of young clusters, without 
accounting for the effects of dynamical evolution. 



The global baryon efficiency can be related to the com- 
monly used specific frequency, S N = iV GC 10°- 4(Mv+15) . 
Taking the mean cluster mass, 2 x 10 5 M Q , and assum- 
ing a mass-to-light ratio for old clusters, M/Ly = 3, we 
obtain Sn = 1.2 x 10 3 M GC /M*. The stellar galaxy mass 
M* is not a very good proxy for the baryon mass Mb at 
high redshifts, but it can be used for the comparison at low 
rcdshifts. Thus, our efficiency £q C s=s 3 x 10~ 4 corresponds 
to Sn ~ 0.4, which is indeed observed for Sc galaxies. 

It is interesting also that the mass of the globular clus- 
ter population and the maximum cluster mass in a given 
region strongly correlate with the local average star for- 
mation rate density: M max oc Egpp ±0 ' 07 and Mqc oc 
Egpp ±0 ' 06 at z = 3.3, where the masses and star for- 
mation rate were estimated taking into account clusters 
with M > 5 x 10 4 M Q and stellar particles younger than 
5 x 10 7 yr and averaging over the cells of 7.7 physical kpc. 
Each averaging cell therefore represents a different progen- 
itor galaxy in the simulation. A similar correlation was 
reported for the observed nearby galaxies (Larsen 2002). 
The interpretation of this correlation is straightforward in 
our model. The star formation rate depends sensitively 
on the mass fraction of gas in cold high-density star form- 
ing regions (Kravtsov 2003b). The massive clusters in our 
model are also assumed to form in such regions. Thus, 
both the star formation rate and the mass of the globular 
cluster population are controlled by the amount of gas in 
the densest regions of the ISM. 

3.4. The mass, size, and metallicity distributions of the 
model clusters 
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Fig. 6. — The mass function of globular clusters formed within 
the parent halos of mass < 3 X 10 10 M© (dotted histogram) and 
> 3 X 10 10 M (solid histogram) at z = 3.3. 



Figure 5 shows the mass function of the model clusters 
identified in all output epochs prior to a given redshift. At 
all epochs the distribution is well described by a power- 
law dN/dM oc M~ a at M > 10 5 M Q . The slope a evolves 
slowly and saturates at a = 2.05±0.07 for z < 4. Note that 
although the physical resolution of the simulation changes 
somewhat with time, the subgrid prescription based on the 
physical threshold density ensures that the cluster prop- 
erties are not affected. In §6.3 we show that the cluster 
mass is a fraction of the mass of the central cell of the 
parent molecular cloud that depends only on the physical 
density of the gas. The mesh in our simulation is refined 
in a quasi-lagrangian fashion, so as to keep the same gas 
mass within a cell, and thus the cell masses and cluster 
masses are always similar. 

In this and subsequent figures we discard the most mas- 
sive cluster identified at the center of the most massive 
disk. Such cluster would be identified in observations as 
a compact galactic nucleus rather than a distinct globular 
cluster. 

The simplicity of the power-law shape of the mass func- 
tion is deceiving, as the mass functions of clusters in in- 
dividual galaxies exhibit a variety of shapes. Figure 6 
shows, for example, that small halos form only clusters 
with M < 10 6 M Q with the mass function slope steeper 
than a = 2, while large halos form more massive clus- 
ters with a shallower mass function. Remarkably, the con- 
volution of the halo distribution with the distribution of 
clusters within each halo produces the seemingly invariant 
power- law mass function of globular clusters. We investi- 
gate the origin of the mass function in detail in § 4. 

In contrast to this power law function, the mass function 
of old Galactic and most extragalactic globular clusters has 
traditionally been described by a bell-shape (Gaussian) 
function. The mass-to-light ratio for old stars is approxi- 
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Table 1 

The quartiles of the size and metallicity distributions of the model and Galactic globular clusters 



R h {25%) R h (50%) R h (75%) [Fe/H](25%) [Fe/H](50%) [Fe/H](75%) 



z = 10 1.9 2.0 2.3 -2.7 -2.6 -2.5 

2 = 7 2.0 2.4 2.7 -2.2 -1.8 -1.6 

2 = 4 2.0 2.4 2.9 -1.8 -1.5 -1.2 

2 = 3.3 2.0 2.3 2.9 -1.7 -1.4 -1.1 

MW, [Fe/H]<-1 1.2 2.8 4.8 -1.8 -1.6 -1.4 

MW, all 0.7 2.4 4.0 -1.7 -1.4 -0.7 



mately constant and it is appropriate to use the luminosity 
function as a proxy to the mass function of globular clus- 
ters. Our results can be reconciled with observations after 
taking into account the effects of subsequent dynamical 
evolution of the model clusters. Sophisticated models of 
the dynamical evolution, started by Spitzer and collabora- 
tors in the 1970s (see Spitzer 1987) and refined in the 1990s 
(e.g., Chernoff & Weinberg 1990; Gnedin & Ostriker 1997; 
Murali & Weinberg 1997; Vesperini & Heggie 1997; Gnedin 
et al. 1999), have shown that tidally-truncated clusters un- 
dergo secular mass loss. The main processes shaping the 
GCMF are the evaporation through two-body relaxation, 
stellar mass loss, tidal shocking by the host galaxy, and 
stronger tidal truncation due to dynamical friction. Fall 
& Zhang (2001) have demonstrated that these processes 
naturally transform the initial power law into the observed 
truncated mass function of old globular clusters. Using a 
simple application of the above results to our model clus- 
ters, we estimate that the dynamical evolution will pro- 
duce the peak and dispersion of the mass function of sur- 
viving clusters in agreement with observations. A more 
detailed analysis, including the orbits of clusters in merg- 
ing galaxies, will be done in a subsequent study. 

Figure 7 shows the distribution of the half-mass radii, 
calculated according to our subgrid model (eq. [4]). As 
are the masses, the cluster sizes do not vary systematically 
with the formation redshift and have consistent distribu- 
tions at all epochs. In this and the following figure we 
plot only the massive clusters, M > 10 5 M Q , expected to 
survive the dynamical evolution. 

The model size distribution is generally similar to the 
observed sizes of the Galactic globular clusters (c.f., van 
den Bergh 1996). The differences at the smallest and 
largest ends can be due to the following effects. The adia- 
batic expansion condition may not apply to the most mas- 
sive clusters which would then be larger. The dynamical 
evolution effects, on the other hand, would shrink the clus- 
ters and fill the range R < 1 pc. Also, some of the sur- 
viving clusters with M < 10 5 M may contribute to the 
smallest size bin as well. 

By construction, our constant density subgrid criterion 
leads to the correlation between cluster masses and sizes, 
R oc M 1 / 3 . Recent observations (Zepf et al. 1999; Larsen 
2004) suggest that the sizes of young massive star clusters 
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Fig. 7. — The size distributions of globular clusters at four suc- 
cessive epochs, using only the massive clusters, M > 10 s Mq . 



are almost independent of their masses and show weaker 
correlation, R oc M 01 . This implies either that the mas- 
sive star clusters form with intrinsically higher densities or 
that the low-mass clusters expand more than we assumed 
following the loss of the remaining gas in the star form- 
ing cores. For example, Ashman & Zepf (2001) suggested 
that the formation efficiency e increases with the binding 
energy of the molecular cloud. It might therefore be useful 
to define two formation efficiencies: the mass conversion 
efficiency sm in equation (3) and the expansion efficiency 
€r in equation (4). As we have argued in §2.4, cm cannot 
be much different from 0.5 — 0.6, but €r can, in principle, 
scale with the cluster mass. The observed trend can be 
explained by oc M 2 . We do not discuss these compli- 
cations further in this paper but they should be addressed 
in future work. 

Figure 8 shows the distribution of cluster metallicities, 
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Fig. 8. — The mctallicity distributions of globular clusters at four 
successive epochs, using only the massive clusters, M > 10 5 Mq . 



which is remarkably similar to the metal-poor part of the 
Galactic cluster distribution (see Table 1). Interestingly, 
we do not find any correlation between the mass and metal- 
licity of globular clusters. The Spearman rank correlation 
coefficient for the cumulative distribution at z w 3.3 is 
r s = —0.06, which is consistent with no correlation. Simi- 
larly, there is no correlation for the Galactic GCs (Djorgov- 
ski & Meylan 1994). The absence of the mass-metallicity 
correlation may be due to a wide range of gas metallici- 
ties in the star forming regions. As we mentioned above, 
this also explains the lack of a well-defined age-metallicity 
correlation. 

At the high redshifts considered here, most of the metals 
are contributed by SNII which underproduce iron com- 
pared to the a-peak elements. In order to calculate the 
fraction of ejected metals contributed by iron, ?7Fe, we use 
the iron yields of Woosley & Weaver (1995), integrated 
over the Miller & Scalo (1979) IMF in the range 12 to 40 
Mq. For the intermediate 'case B' yield models and the as- 
sumed solar ratio of iron to hydrogen of 1.8 x 10 -3 by mass 
(Anders & Grevesse 1989), we obtain t]f c w 0.6. This es- 
timate is consistent with the enhanced ratios [a/Fe] « 0.3 
observed for the globular clusters in the Galaxy (Carney 
1996; Lee & Carney 2002; Smith et al. 2002), M87 (Cohen 
et al. 1998), and M49 (Cohen et al. 2003). The conversion 
from the total metallicity due to SNII to the iron abun- 
dance is [Fe/H] = \g(-q Fe Z n /Z Q ). 

Table 1 compares the medians and the 25% and 75% 
quartiles of the size and metallicity distributions of the 
model clusters against the corresponding distributions for 
the Galactic GC (Harris 1996). There is good agreement 
between our predictions and the observations. It is plau- 
sible that small discrepancies for the smallest sizes and 
highest metallicities can be explained by subsequent evo- 
lution of the cluster population at z < 3, which we discuss 
in §6.1. Note that the evolution can also modify the quar- 



tiles of the distribution. We therefore show in Table 1 both 
the quartiles for all Galactic clusters and for the clusters 
with metallicities in the same range as in our model (i.e., 
[Fe/H] < — 1]). It is also interesting that in agreement 
with observations no clusters are formed with very low 
(Pop III) metallicities. 

Note that at all epochs the dynamical time of the parent 
molecular cores is very short (~ 10 6 yr), which means that 
the galactic gas is pre-enriched even before the first clusters 
form. It follows then that the oldest globular clusters do 
not contain the oldest stars in the Galaxy. 

The lack of the metal-rich clusters in our model com- 
pared to the Galactic clusters (c.f. Table 1) is likely to be 
explained by the clusters forming in the higher-metallicity 
gas at z < 3 (see § 6.1), not accounted for in our analysis. 
It should be noted, however, that a significant spread in 
metallicity distributions exists for different galaxies (e.g., 
Harris 2001) and certain differences with the simulated 
system are expected. This issue can be addressed in fu- 
ture studies by simulating globular cluster systems in a 
number of galaxies. 

4. THE ORIGIN AND UNIVERSALITY OF THE GLOBULAR 
CLUSTER MASS FUNCTION 

One of the most important characteristics of globular 
cluster systems is the mass function (GCMF). Interest- 
ingly, the mass function derived in the simulation is sim- 
ilar to the mass function of molecular clouds in our and 
external galaxies, dN/dM oc M~ a with a = 1.4 — 2 (e.g., 
Solomon et al. 1987; Wilson et al. 2003). It is also similar 
to the high-redshift mass function of dark matter halos in 
the hierarchical CDM cosmology (Press & Schechter 1974; 
Sheth & Tormen 1999). In this section we investigate the 
origin of GCMF in relation to the mass function of giant 
molecular clouds and parent halos. 

Gncdin (2003) used a simple semi-analytic model in 
which a single massive cluster dominates the mass of the 
globular cluster system within a progenitor halo: M max < 
M GC « SQ C f h M h (c.f. § 3.3). The model implies that the 
shape of the high-mass tail of the GCMF simply reflects 
the shape of the mass function of progenitor halos of the 
Milky Way at high redshifts. This direct connection be- 
tween the cluster and halo mass functions is due to the 
assumption that GC properties are determined solely by 
the mass of their parent halo. This key assumption can be 
tested against the results of our simulation. 

We find indeed that the most massive cluster contributes 
a significant fraction of the total cluster mass. The average 
for all halos is M max « 0.6 Mqc- In § 3.3 we have shown 
that Mqq is roughly proportional to the parent galaxy 
mass, and therefore a similar relation exists for the most 
massive clusters: 

/ M \ 1-29±0.12 

M max = 2.9xlO 6 M (^ T ^-j . (8) 

However, a significant scatter around this average relation 
is such that for a given halo the masses of individual clus- 
ters can vary by a factor of three. Thus, the mass function 
of globular clusters does not follow directly from the mass 
function of their parent halos. 

Instead, we find that the overall shape of GCMF in our 
model is determined by both the mass function of progen- 
itor halos and the mass function of molecular cloud cores 
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Fig. 9. — Probability distribution functions of the gas density 
(fraction of the total volume occupied by the cells in a given density 
range) for the highest level of refinement at z = 7, 4, and 3.3. The 
dotted lines show the log-normal fits to the highest-density tails of 
the distribution, in the range 1 — 10 3 Mq pc -3 . The straight lines 
show the slopes dV/dlgp oc p _1 and p -15 . 



within individual halos. The latter, in turn, is determined 
by the structure of the galaxy disks, and in particular by 
the density probability distribution function (PDF). 

The shape of the PDF has been studied in several nu- 
merical simulations of the turbulent ISM. A log-normal 
distribution is thought to be a generic feature of isother- 
mal turbulent flows (e.g., Vazquez-Semadeni 1994; Padoan 
et al. 1997; Vazquez-Semadeni et al. 2000). For non-isothermal 
supersonic flows Scalo et al. (1998) found that the PDF is a 
power law of density, although Nordlund & Padoan (1999) 
argued that even in this case the PDF can be described 
equally well by a power law or a log-normal function. The 
log-normal shape is likely to be due to the chaotic nature 
of the supersonic turbulent flows, characterized by numer- 
ous random convergent flows and shocks. The evolution 
of individual gas elements can be thought of as a random 
walk in density leading to the log-normal equilibrium dis- 
tribution (Elmegreen 2002). 

For the heating/cooling rates adopted in our simulation 
the gas at densities p g > 5 M Q pc" 
the lowest allowed temperature (T„ 
therefore nearly isothermal. Accordingly, we expect the 
PDF in the simulation to be log-normal. The width of the 
log-normal PDF is proportional to the rms Mach number 
of gas clouds (Padoan et al. 1997), which increases with 
time in the hierarchically assembled galaxies. Thus we 
expect the PDF to widen with decreasing redshift. 

Figure 9 shows the density PDF measured in the simu- 
lation at z = 7, 4, and 3.3, using only the highest refine- 
ment level (I = 9) cells which cover the galactic disks and 
include the sites of GC formation. The log-normal distri- 
bution provides a very good fit to the high density tail of 
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PDF at p > 1 M pc~ 
Kravtsov 2003b): 
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(see also Wada & Norman 2001; 
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The characteristic density po and the dispersion a p of the 
log-normal PDF vary with redshift. The characteristic 
density decreases from lgpo ~ 1-8 at z = 8 to lgpo ~ 0.08 
at z — 3.3, while the width of the distribution a p increases 
from 0.46 to 0.85 at the same redshifts. The evolution in 
this redshift interval can be fit by lgpo = 3.30 — (1.41 ± 
0.02) [10/(1 + ;?)] and a p = 1.23 - (0.83 ± 0.05)[(1 + z)/10]. 

Figure 9 shows also that over a limited density range, 
1 < lg p < 3, the PDF can be described by a power-law 

cx p- (10) 

d\gp 

with n = 1.08 ± 0.06 (for z = 4). Over a wider range of 
densities, however, the log-normal function is a somewhat 
better description of the PDF. For instance, the power-law 
slope becomes steeper with increasing density: we find 
n = 1.26 ± 0.08 for Igp > 1.5 and n = 1.41 ± 0.13 for 
\gp > 2. The latter range includes the molecular cloud 
cores in which the massive clusters (M > 10 5 M ) form in 
our model. 

We should note that the molecular clouds in the sim- 
ulation arc only marginally resolved. Their temperature 
and structure depend on the uncertain physics such as the 
presence of dust, cooling due to metals, UV heating by 
nearby stars and radiative transfer. The form of the PDF 
can depend on these properties and at this point we can- 
not reliably differentiate between log-normal and power 
law mass functions. However, the difference between the 
two at the interesting densities is quite small and is not 
critical for our discussion. 

Globular clusters in our model form only in the highest- 
density cores of the identified molecular clouds which rep- 
resent a subset of all high-density cells. Nevertheless, we 
find that the PDF of the cores (or density peaks) is simi- 
lar to the overall cell PDF at lgp > 0. The Kolmogorov- 
Smirnov test of the unbinncd probability distributions shows 
that the differences between the core PDF and the cell 
PDF arc not statistically significant. At all epochs the 
probability that the two PDFs are drawn from the same 
distribution is at least 13%. This indicates that the den- 
sity PDF of the cores is also described by the same log- 
normal distribution. The density of each core determines 
the mass of the globular cluster it hosts. The density PDF 
of the cores thus determines the mass function of globular 
clusters. 

Given our subgrid model, the cluster mass scales with 
the core gas density as M oc p 3 / 2 (eqs. [3] and [4]). For a 
power-law density PDF the expected cluster distribution 
is dN/dM oc M- 1 - 2 "/ 3 , or a = 1 + 2n/3. For n = 1 this 
gives the slope a = 5/3 « 1.7, while for n = 1.4 ± 0.13, 
appropriate for the highest-density tails of the PDF and 
the massive clusters, a = 1.94 ± 0.09, in good agreement 
with the mass function slope seen in Figure 5. Therefore, 
we can expect a relatively shallow, a « 1.7, mass function 
for the small-mass clusters forming in lower-density cores 
and the steep, a « 2, mass function for massive globular 
clusters forming in the densest regions of the disk. The 
range of slopes derived in our model is in close agreement 
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Fig. 10. — Mass function of globular clusters forming at two 
epochs (histograms) with the superimposed fit (smooth solid and 
dashed lines) from the density PDF (eq. [11])- The error bars 
represent Poisson errors oc N 1 / 2 . The dashed line shows the mass 
function calculated from the density PDF of the lower, I = 8, level 
of refinement normalized to the same total number of clusters. 



Fig. 11. — The mass function of young (age 25 < t < 160 Myr) 
star clusters (solid circles with error bars) in NGC 4038/9 (the "An- 
tennae"), as derived by Zhang & Fall (1999), and the best fits of 
the power law (dotted line) and log-normal (solid curve) mass dis- 
tributions. 



with observations (Elmegreen & Efremov 1997; Zhang & 
Fall 1999; de Grijs et al. 2003; Anders et al. 2004). 

From our analysis above it is clear that the power-law 
is not a unique description of the mass function. The log- 
normal shape of the density PDF implies the log-normal 
GCMF: 
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with 



cm = 2 CT p 



and M = M(pq) using equation (3). Fig- 
ure 10 shows the mass function of clusters at z = 8 and 
z = 4 along with the log-normal function calculated from 
the best fit to the density PDF. In other words, the smooth 
lines in the figure represent not fits to the mass function, 
but the fits to the PDF converted to the mass function us- 
ing our subgrid model. The log-normal mass function may 
thus be a reasonable choice in fitting observed luminosity 
and mass functions of young clusters when a simple power 
law does not provide a good fit. 

The parameters of the distribution, M and <7m , follow 
directly from the parameters of the density PDF. As red- 
shift decreases, the characteristic peak p decreases and 
the dispersion a p increases. While their exact values at a 
given redshift may be specific to our simulated galaxy, i.e. 
depend on the environment, the anti-correlation between 
the parameters may be more general. We find the follow- 
ing relation which can be tested by future simulations and 
observations: 

lgM = (6.2 ± 0.5) - (2.8 ± 0.5)cr M . (12) 
These results apply to the mass function of young stellar 
clusters not modified by dynamical evolution. 

In order to illustrate that the log-normal distribution 
can fit the observations, we plot on Figure 11 the mass 



function of young (age 25 < t < 160 Myr) stellar clusters 
in NGC 4038/9 (the "Antennae"), as derived by Zhang & 
Fall (1999). This figure shows that the log-normal function 
describes the observed mass function as well as the power 
law. The best fit power-law slope is a = 2.01 ± 0.07 with 
the reduced % 2 = 1.17, while for the log-normal fit \ 2 = 
0.68. The parameters of the fit, lgM = 3.7 ± 0.4 and 
<7m = 0.73 ± 0.12, are in reasonable agreement with the 
derived relation (12), especially if we take into account the 
difference in rcdshifts, metallicity, and environment. 

The results presented in this section indicate that the 
observed luminosity and mass functions of young clusters 
could be described by a power-law due to the limited range 
of luminosities typically probed in observations (^2 — 4 
magnitudes or only about a factor of 10 — 40 in mass: 
Elmegreen & Efremov 1997; Whitmore 2000). If our re- 
sults are correct, the prediction would be that the power- 
law slope a should steepen at the highest cluster masses: 
M > 1O 7 M . Interestingly, the lo g-normal mass function 
implies that there exists a maximum cluster mass at any 
given epoch. This is an important conclusion as it may 
explain the characteristic masses of clusters in the Local 
Group. 

5. NUMERICAL CONVERGENCE 

The resolution of the simulation in any numerical study 
invariably places constraints on the validity of the results. 
The spatial and mass resolution of the gas element in the 
adaptive mesh simulation is set by the maximum level of 
refinement. In our simulation l max = 9 was determined 
by the available computational resources. We can check, 
however, how the results would change if we limited the 
maximum level to I = 8. This corresponds to factor of 8 
lower mass resolution and a factor of two larger cells. 
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We find that the density PDF of the I = 8 cells can be 
well fit by a log-normal function but with a correspond- 
ingly lower characteristic density p n . The dispersion a p 
is consistent with that for I = 9 level within the errors. 
The mass function resulting from the / = 8 density PDF 
is therefore somewhat steeper than for the / = 9 cells, be- 
cause M is lower and the same mass interval of globular 
clusters falls on the steeper part of the log-normal func- 
tion. The difference, however, is only apparent at very low 
masses (M < 3 x 10 4 M Q ). 

The dashed line in Figure 10 shows the cluster mass 
function expected for the density PDF of the I = 8 cells 
normalized to the same number of clusters as the I = 9 
mass function. It fits the histogram of the model clusters 
equally well. The deviation only becomes significant at 
M < 3 x 10 4 M Q . Thus, the shape of the derived mass 
function converged at masses higher than 3 x 1O 4 M . 

Our assumption that each molecular cloud forms only 
one cluster may also affect the low mass end of the mass 
function. If in reality the cloud fragments into several 
self-gravitating cores, then smaller clusters may form on 
the periphery of the cloud in addition to the larger central 
cluster. However, as we show in the next section, these ad- 
ditional small clusters are likely to be quickly dissolved, so 
that in the end the mass function of the surviving massive 
clusters remains the same. 

6. DISCUSSION 

6.1. Evolution of Globular Clusters at lower redshifts 

In the preceding sections we analyze the formation of 
globular clusters at z > 3. Although the limited computa- 
tional resources did not allow us to continue the simulation 
at the same resolution to lower redshifts, in this section we 
conjecture on the possible evolution of cluster population 
at later epochs. 

Analysis of our simulations hints that massive molecular 
clouds needed for cluster formation are built in gaseous 
spiral arms. Their formation may be enhanced in mergers 
between gas rich gas disks. For example, as we note in 
the caption of Figure 1, the gas disk shown in this figure 
is actually in a state of very active accretion and merging. 
For instance, the two cold dense clumps clearly seen in the 
temperature map are two satellite galaxies in the process 
of merging with the disk. The mass of the parent halo of 
the disk in Figure 1 increases by a factor of twenty between 
redshifts of 10 and 4 (a period of only ~ 1 Gyr! See § 2.2). 

It is commonly thought that galaxy mergers create con- 
ditions conducive to bursts of star formation and star clus- 
ter formation, as evidenced by the starbursting galaxies 
(Kennicutt 1998; Larson 2002). In particular, mergers 
stir and compress the interstellar gas creating the high- 
pressure environments in which dense clouds and massive 
stellar clusters can form (Elmegreen 2002). Without the 
mergers, star formation proceeds in a quiescent mode (e.g., 
Abadi et al. 2002). 

The high rate of accretion and merging cannot be main- 
tained at lower redshifts. Statistically, the CDM halo 
merger rate at z < 4 decreases rapidly as (1 + z)~ 2 - 5 
(Gottlober et al. 2001). In addition, as galaxies evolve, 
most of the gas will be converted to stars so that the 
large reservoir of gas needed to build up giant molecular 
clouds may not be available. Therefore, if galaxy mergers 
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Fig. 12. — The number of clusters (triangles) and the total mass in 
clusters (circles) formed at each simulation output according to our 
model. The number of sites of cluster formation reaches a plateau 
at z ai 4 and may decline at lower redshifts. 



are connected to globular cluster formation, a high merger 
rate between gas rich galaxies at z > 3 would lead to an 
almost continuous cluster formation. At lower redshifts, 
on the other hand, mergers become rare and the merging 
galaxies are gas deficient compared to their high-redshift 
progenitors. Most of the subsequent formation of stellar 
clusters may thus be limited to a single last major merger 
event. This, for example, could explain the bimodality of 
cluster colors observed in many elliptical galaxies. 

Finally, somewhat paradoxically, the globular clusters 
formed at high redshifts may have a significantly higher 
chance of survival until the present than clusters formed at 
later epochs. The clusters in our model form in extremely 
high-density environments within the galactic disks. Strong 
tidal forces in such regions are likely to disrupt clusters 
quickly, unless they are located at the very center of the 
parent galaxy or are ejected from the disk by a dynamical 
process shortly after formation. As only a fraction of clus- 
ters form in the centers of progenitor galaxies, the latter 
mechanism should operate in order for the old metal-poor 
clusters to survive until the present. 

Frequent violent mergers at high redshifts may be just 
such a mechanism. Mergers disrupt the disks of the pro- 
genitor galaxies with their existing stellar populations and 
impart a large amount of orbital energy in the surviving 
clusters. The high energy orbits would allow globular clus- 
ters to spend most of the time in the relatively low-density 
regions of the halo, outside the main disk. Mergers could 
be responsible for a spheroidal distribution of the globu- 
lar cluster systems, even though the clusters actually form 
within the disks of the progenitor systems. At low red- 
shifts, on the other hand, star clusters, even if they con- 
tinue to form, may not be able to escape the disk quickly 
enough to avoid disruption. 
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The effects discussed above may result in the preferred, 
albeit extended, epoch of globular cluster formation at 
zgc ~ 2 — 12 (tec ^0.3 — 3 Gyr for the adopted cosmo- 
logical model). We cannot prove this conjecture directly 
in our present simulations, although one may argue that 
there are indications of this trend at z ~ 3 — 4. Figure 
12 shows that the total mass of the GC population and 
the number of clusters forming at a given epoch increases 
until z « 4 and then saturate at a constant value at lower 
redshifts. 

Note that Figure 12 should be interpreted with caution. 
Since clusters are not formed self-consistently during the 
simulation run time but rather calculated a posteriori, it 
is impossible to compute the rate of their formation. In a 
self-consistent treatment, clusters forming at a given epoch 
may exhaust the supply of gas in the regions conducive to 
cluster formation and prevent the formation of new clus- 
ters in the same place. In our current analysis, on the other 
hand, clusters at subsequent epochs can, in principle, form 
from the same gas. At the high redshifts that we consid- 
ered, the gas distribution changes sufficiently quickly and 
the rate of accretion of new gas is very high, so that the 
epochs separated by tens of millions of years can be safely 
considered independent. This may not be true for nearby 
epochs. The flattening of Mqc an d Nqc in this figure 
can be interpreted as the fact that no new sites for clus- 
ter formation are created. Unfortunately, this is the best 
estimate of the formation rate that we can provide with 
our current simulation. Our future work will address the 
formation of star clusters at low redshifts directly. 

To summarize, we conjecture that the formation of glob- 
ular clusters at z < 2—3 can be associated with the increas- 
ingly rare merger events and would thus be progressively 
episodic. The merger rates derived from cosmological sim- 
ulations (Gottlober et al. 2001) indicate that on average 
galactic halos experience 1 to 4 major mergers with the 
mass ratio > 0.2 at z < 2. Two broad evolutionary sce- 
narios can be envisioned in this picture. If the last major 
merger occurs early (2: m0 rgc > 2), not too far from the 
preferred epoch of GC formation, we would expect a con- 
tinuous change of the cluster properties. If, on the other 
hand, the galaxy experiences the last major merger late 
( Emerge ^ 1), with a substantial gap between z m orgo and 
zqc 7 the distribution of properties of the resulting cluster 
populations may be bimodal. The bi-modality in this case 
would reflect a considerable change in the galactic envi- 
ronment (e.g., the metallicity of gas) during the interval 
Az — zqq — z morgo . This latter scenario is particularly 
relevant to the formation of large elliptical galaxies. 

6.2. Comparison with previous work 

Several recent studies explored the formation of globular 
clusters in the context of hierarchical cosmology using both 
numerical simulations and phenomenological semi-analytic 
models. Here we discuss and compare the specifics and 
results of our study to the previous similar efforts. 

Weil & Pudritz (2001) used a Tree-SPH simulation of 
the rCDM (fl = 1) cosmology to study the large-scale 
distribution of giant gas clouds at z < 1 in a 32h~ 1 
Mpc box. The authors analyzed collisions of baryonic 
clumps (clouds) in small-mass dark matter halos within 
the context of the agglomeration model of Harris & Pu- 



dritz (1994). They found a characteristic power-law spec- 
trum of cloud masses, dN/dM mc cx M" 1 / 7 , similar to the 
mass function of molecular clouds and young star clusters. 
The mass function of clouds and globular clusters in our 
model is consistent with their result. The detailed com- 
parison, however, is not possible as the simulations of Weil 
& Pudritz (2001) did not include star formation and the 
cooling of gas below 10 4 K. In addition, their low spatial 
resolution (~ 1 kpc) prevented any detailed study of the 
inner structure of the clouds and as well as the density dis- 
tribution and mass spectrum of clusters within individual 
galaxies. 

Recently, Bromm & Clarke (2002) used a Tree-SPH sim- 
ulation to study the collapse and fragmentation of gas 
during the evolution a single dwarf galaxy ~ 10 8 M Q at 
z < 24. The simulation assumed that the gas was pre- 
enriched to the metallicity of 10~ 2 Z@, which allowed the 
gas to cool to ~ 1000 K, and used sink particles to follow 
the collapse of the gas in the highest density (n > 10 3 
cm~ 3 ) regions. The authors identified six gas clumps with 
masses in the range 4 x 10 4 — 2 x 10 7 M Q , each associ- 
ated with a separate small-mass dark matter halo. They 
concluded that the characteristic mass of globular clusters 
is determined by the characteristic mass of dark matter 
halos forming at z > 10. 

The implicit assumption behind their conclusion is that 
the conditions for cluster formation exist only at the ear- 
liest stages of galaxy formation, prior to reionization. The 
results of our study show that, although the first globular 
clusters may form at z > 10, the conditions for cluster 
formation become more favorable at lower redshifts when 
more gas accumulates in the disks of the progenitor halos. 
Reionization does not affect significantly the formation of 
clusters in the relatively massive halos (Mh jcj 10 10 M ). 
Again, it is difficult to make a more detailed comparison, 
given the very different setup of numerical simulations and 
physical processes included. Yet, it is worth noting that 
the gas density in the most massive clump in the simula- 
tion of Bromm & Clarke (2002) (see their Fig. 3) is well 
below the density of dark matter. The average gas density 
at the resolution limit in the inner 10 pc is only ~ 10 Mq 
pc , lower than the observed density of globular clusters 
(§ 2.4). The expected mass loss after cluster formation 
would reduce the cluster density even further. It is likely 
that such clumps are still insufficiently dense to form real 
globular clusters. 

In contrast, the gas density at the sites of globular clus- 
ter formation in our simulation is considerably higher (> 
100 M pc -3 ) and is typically at least an order of mag- 
nitude larger than the local density of dark matter. We 
further assume that the clusters form only at densities 
p > 10 4 Mq pc -3 within the collapsing isothermal molec- 
ular cores. The fact that GC formation sites are strongly 
baryon-dominated explains the absence of dark matter ha- 
los around globular clusters. Bromm & Clarke (2002), 
on the other hand, conjecture that dark matter hosts of 
globular clusters dissolve via violent relaxation before the 
cluster forms, while the baryonic cores survive. 

Both of the above studies attribute the shape of the GC 
mass function to the distribution of their parent dark mat- 
ter halos. Our results indicate that the situation is more 
complex. We have shown that the shape of the GCMF is 
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determined both by the mass function of the parent halos 
and the mass distribution of clusters within a single halo 
(§ 4). 

Bcasley et al. (2002) extended the semi-analytical model 
of galaxy formation GALFORM to include a phenomenolog- 
ical prescription for GC formation. The model assumed 
that a constant fraction of stellar mass, e, would be in the 
form of globular clusters. The blue (metal-poor) clusters 
were associated with the quiescent mode of star forma- 
tion, while the red (metal-rich) clusters were assumed to 
form during starbursts. The parameters, e = 0.002 for 
blue clusters and e = 0.007 for red clusters were set to 
agree with the observed color distribution of the elliptical 
galaxy NGC 4472. However, Beasley et al. (2002) trun- 
cated the formation of blue clusters arbitrarily at z = 5 to 
create a distinctly bimodal distribution of colors in their 
model. In addition they find an age-metallicity correlation 
for their model clusters, which is a definite prediction of 
semi-analytical models. In contrast, our simulation follows 
the gas dynamics and metal enrichment self-consistently 
and does not predict any clear age-metallicity relation. 

6.3. Alternative subgrid models 

The results we presented are based on a particular sub- 
grid model. The main features of the adopted model are 
the assumption of the isothermal cloud structure within 
the density peak of each molecular cloud and the assump- 
tion that clusters form with a fixed efficiency at densities 
above a constant density threshold (see § 2.4). 

Here we consider alternative subgrid models and show 
that they can either be reduced to the model we use or 
cannot successfully reproduce the mass function of globu- 
lar clusters. Specifically, we consider the following alter- 
native functions for the cluster formation threshold: (1) 
the Jeans mass, (2) the thermal pressure, and (3) the to- 
tal pressure contributed by thermal and by turbulent mo- 
tions. The size of the star forming core, and therefore the 
cluster mass, is determined differently by a threshold value 
of each of these functions. 

In general, an alternative subgrid model can have a 
threshold parameter that varies with the radius differently 
from the isothermal gas density profile, p g oc r~ 2 . It is 
likely, however, than in an isothermal cloud any function 
of interest would be a power law, /(r) oc r~ a . The star 
formation radius corresponding to the threshold value / cs f 
is then determined by f(R cs f) = /csf- 

Since in the simulation we can only measure the parame- 
ter / ce n averaged over the cell volume, we need to consider 
the volume average 

/av(< r) = ^ f f(r)Avr 2 dr = ^-/(r). (13) 
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Fig. 13. — The build-up of the initial mass function of globular 
clusters based on the turbulent pressure criterion. Dotted, dashed, 
long-dashed, and solid histograms show cumulative distributions at 
z = 8, 7, 4, and 3.3, respectively. The straight dashed lines show two 
power-laws, TV oc M~ Q . Note that the mass function is significantly 
shallower than that with the density criterion, a = 2. 



where M ce n = (47r/3)p ce iiP 3 ell is the amount of gas within 
a sphere embedded into the cell. Including the efficiency 
of star formation e, the stellar mass is = eM(R s {) and 
stellar radius is P* = P s f/e. 

The Jeans mass is a critical mass of a cloud of density 
p g and temperature T that is unstable to gravitational 



For isothermal molecular 
, a threshold in Mj 



instability, Mj cx T 3 / 2 p g 1/2 . 
clouds with T w const and p g oc 
reduces to a corresponding threshold in p g , i.e. our original 
model. Similarly, the thermal pressure criterion P tn oc p g T 
reduces to the density criterion. 

The only variable independent of the gas density, at least 
in principle, is the turbulent velocity dispersion, a. Turbu- 
lence in the galactic disks is created by gravitational mo- 
tions on scales larger than a single cell and reflects large- 
scale flows around the cloud. In the simulation, the turbu- 
lent pressure, Pturb = Pg ' 2 , typically dominates over the 
thermal by one or two orders of magnitude, so that the to- 
tal pressure P « -Pturb- The velocity dispersions and sizes 
of the Galactic giant molecular clouds satisfy the follow- 
ing relation: a oc r 1 / 2 (Larson 1981; Brunt & Heyer 2002). 
If the turbulent velocity dispersion scales with the radius 
inside the cell as a oc r 1 / 2 , while p g oc r~ 2 , then Par -1 . 

A reasonable fiducial value for the threshold pressure, 
P cs {, can be obtained from observations of the Galactic 
globular clusters. The median value of the pressure at the 
half-mass radius is P b s ~ 10 9 K cm~ 3 . The observable 
pressure in our model can be derived taking into account 
the post-formation expansion of the cluster. Analogously 
to the observable density (eq. [5]), we find P = |e 4 P cs f. 
Setting P = P bs, w e obtain the star formation threshold 
P csf = 5 x 10 9 K cm" 3 = 10 6 Mq pc~ 3 km 2 s~ 2 . 
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The velocity dispersion in a cloud is calculated directly 
from the gas velocities in the simulation. First, the cloud 
is centered on the cell of highest density. Then, the center- 
of-mass motion, the net radial motion, and the solid-body 
rotation velocity are subtracted from the velocities of the 
cloud cells. The remaining velocity field is free of global or- 
ganized motion. The resulting velocity dispersion is close 
to isotropic, as expected for a clean turbulent field. We 
use these turbulent dispersions, a, to calculate the tur- 
bulent pressure in the central cloud cell. Then we apply 
equations (15) and (16) to calculate the predicted cluster 
sizes and masses, using / ceU = p C eii^ ell , / cs f = P cs f, and 
a = 1. 

Figure 13 shows the cluster mass function calculated for 
this subgrid model, which should be compared to Figure 
5. It is clear that the turbulent pressure criterion predicts 
a much shallower mass function than that in the density 
threshold model. The power law slope is a < 1.5. Such 
shallow slope is inconsistent with observations. 

The value of the slope can be easily understood. Accord- 
ing to equation (16), the cluster mass is a certain fraction 
of the cell mass. This fraction depends on (Pceii/Pcsf) 1 ^ 11 , 
and in this case a = 1. If the fraction can be expressed as 
some power of the cell mass, M cc n, then the cluster mass 
function can be derived from the mass function of cen- 
tral cloud cells (uniquely determined in the simulation), 
and the transformation would depend on the particular 
subgrid model. Most of the densest cells are at the last 
refinement level and all have the same volume at any given 
epoch, so that M cc \\ oc p cc n. In agreement with the den- 
sity PDF (eq. [10]), the distribution of the central cells is 
roughly described by a power law 
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for Mcdi > 10 7 M , which are the cells harboring most of 
the massive clusters. 

Let us now derive the transformation from the cell mass 
to the cluster mass. The turbulent velocity does not in 
general scale with the cell mass, since it is caused by large- 
scale flows, but we still find a fairly robust correlation 
in the expected sense, a 2 oc M cc \\, to within 0.3 dex in 
lgc 2 . Therefore, M* oc Pccii-Pccii oc M^ ell . Substituting 
this relation to equation (17), we find 



dN 

dM* 



oc Ml 



3/2 



(18) 



This is comparable to the slope seen Figure 13 and is sig- 
nificantly shallower than the observed slope. 

The general formalism of calculating the cluster mass 
function developed here gives us a powerful tool to eval- 
uate the alternative subgrid models. Unless the transfor- 
mation from the cell mass to cluster mass is the same as 
for the density criterion (M* oc leading to a w 2), 

the predicted cluster mass function would deviate from the 
observed. The criterion based on the turbulent pressure, 
at least in principle independent of the density criterion, 
fails to reproduce the observed mass function of young star 
clusters. 

Of course, our prescription based on a simple threshold 
value of a spherically-symmetric function does not capture 
all the details of the molecular cloud physics, such as the 
filamentary internal structure, dust, and magnetic fields. 



These complications are beyond current models of galaxy 
formation. However, within the scope of our spherically- 
symmetric approach to the molecular clouds, the cluster 
formation criterion based on the density threshold seems 
to reproduce the observations best. 

7. SUMMARY 

We have presented a study of globular cluster formation 
at z > 3 using a very high-resolution cosmological simula- 
tion. The clusters in our model form in the high-density 
isothermal cores of giant molecular clouds in dense gaseous 
disks of high-redshift galaxies. The properties of globular 
clusters are estimated using a simple physically-motivated 
subgrid model. Many of the observed properties of globu- 
lar clusters are reproduced naturally, without fine-tuning. 
Our main conclusions can be summarized as follows. 

• Our results indicate that the first globular clusters 
in the Galaxy should have formed around z « 12 
and cluster formation continued at least until z w 3. 

• Most globular clusters in our model form in halos of 
mass > 10 9 Mq. The spatial distribution of these 
halos at z > 3 is highly clustered (biased) with re- 
spect to the overall distribution of matter (§ 3.1). 
The high spatial bias of their parent halos explains 
the present more concentrated radial distribution of 
globular clusters relative to dark matter. 

• Within the progenitor systems, globular clusters 
form in the highest-density regions of the disk: the 
cores of molecular clouds. In the most massive disk 
in the simulation, the newly formed clusters trace 
the spiral arms and the nucleus similarly to the 
young stellar clusters observed in merging and star- 
bursting galaxies. 

• The mass function of globular clusters at birth can 
be approximated by a power-law dN/dM oc M~ a 
with a w 2, in good agreement with observations of 
young star clusters in the Antennae. It can also be 
well approximated by the log-normal function. The 
shape of the mass function is determined both by 
the mass function of parent galaxies and the mass 
distribution of molecular cloud cores within each 
halo. Although the physical processes governing the 
evolution of the halos and molecular clouds are com- 
pletely different, their mass distributions can both 
be approximated by the power-law functions with 
aw 1.5 -2. 

• The halo mass function arises during the hierar- 
chical build-up of structures in the expanding uni- 
verse from the initial random perturbations. The 
mass function of the molecular clouds cores, on the 
other hand, samples the underlying density proba- 
bility distribution function (PDF) of the gas in the 
galactic disks. 

• The local efficiency of globular cluster formation 
within the parent molecular cloud is M/M mc w 
10~ 3 , with considerable scatter. The global effi- 
ciency, the ratio of the total globular cluster mass to 
the baryonic mass of the parent galaxy, is Mqc /M^ r 
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(2 — 3) x 10~ 4 . The mass in globular clusters scales 
with the total galaxy mass as Mqc oc M^ 1 (eq. 
[6]). These values for the formation efficiencies are 
in general agreement with observations (Harris & 
Pudritz 1994; McLaughlin 1999). 

• We find that the mass of the globular cluster pop- 
ulation and the maximum cluster mass in a given 
region strongly correlate with the local average star 
formation rate density: M max oc Egpp and 
M GC oc £°j^ ±0 - 06 at z = 3.3. A similar correla- 
tion exists for the observed nearby galaxies (Larsen 
2002). The correlation arises because both the star 
formation rate and mass of the globular cluster pop- 
ulation arc controlled by the amount of gas in the 
densest regions of the ISM. However, it is not clear 
whether this correlation is general or applies only 
to gas rich starbursting environments. 

• Our model predicts the lack of clear age-metallicity 
and mass-metallicity correlations, at least for the 
clusters with [Fe/H] < —1. Although there is an 
overall trend of increasing the average metallicity 
with time, a significant spread of metallicities exists 
among different progenitor galaxies and within the 
interstellar medium of each galaxy. The field stars 
and stellar clusters forming at a given epoch in the 
simulation exhibit scatter in metallicity of up to two 
orders of magnitude. 

• The distribution of sizes and metallicities of the 
massive (M > 10 5 M Q ) globular clusters match 
those of the Galactic globulars, with the exception 
of the largest size and the highest metallicity tail. 
It is plausible that these discrepancies can be rec- 
tified by dynamical evolution and the continuing 
formation of globular clusters at z < 3. 

The simulation results presented in this paper, directly 
address only formation of clusters at z > 3. We conjecture 
(§ 6.1) that the overall efficiency of globular cluster forma- 
tion and chances for their survival are significantly reduced 
at lower redshifts with most clusters at z < 2 forming in 
rare gas-rich galactic mergers. Although we cannot prove 
this with our current simulations, we argue that the pre- 
ferred formation epoch of globular clusters surviving until 
the present is z ^ 3 — 5 (or t ^ 1 — 2 Gyr in the adopted 
cosmology) when the gas supply is abundant in the disks 
of the progenitor halos and the merger rate of the progen- 
itors is high. All old globular clusters thus would appear 
to have similar ages. 

Although the high-redshift globular clusters form in dense 
gaseous disks, most the high-z disks arc expected to merge 
with the Milky Way and disrupt by z = 0. The disrupted 
systems form diffuse dark matter halo and contribute to 
the stellar halo of the host. Their clusters would share the 
fate of the stripped stars of the disrupted galaxies that 
build up galactic stellar halo and should therefore have 
present-day spatial distribution similar to that of the stel- 
lar halo. 

In future work, we plan to study the details of the dy- 
namical evolution of the globular cluster population by 
extending our current simulation to z — 0. We will incor- 



porate our subgrid model of star cluster formation directly 
in the high-resolution simulations of galaxy formation. 

Our present results are encouraging and demonstrate 
that globular clusters with properties similar to the ob- 
served clusters form naturally within young z > 3 galaxy 
disks in the standard ACDM cosmology. The conditions 
for cluster formation appear to be widespread at z ~ 5, 
coinciding with the peak of global star formation rate. As 
the formation of first stars marks the end of cosmic dark 
ages (Rees 1997), the formation of globular clusters marks 
a veritable stellar renaissance of the Universe. 
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